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Abstract 

At very high energies, in the infinite momentum frame and in light cone gauge, a hard scale 
proportional to the high parton density arises in QCD. In an effective theory of QCD at 
small x, this scale is of order ag/i, where fi is simply related to the gluon density at higher 
rapidities. The ah initio real time evolution of small x modes in a nuclear collision can be 
described consistently in the classical effective theory and various features of interest can be 
studied non-perturbatively. In this paper, we discuss results from a real time SU(2) lattice 
computation of the production of gluon jets at very high energies. At very large transverse 
momenta, kt > fi, our results match the predictions from pQCD based mini-jet calculations. 
Novel non-perturbative behaviour of the small x modes is seen at smaller momenta k t ~ aj/j. 
Gauge invariant energy-energy correlators are used to estimate energy distributions evolving 
in proper time. 



1 Introduction 



It is of considerable theoretical and experimental interest to understand the collisions 
of nuclei at ultrarelativistic energies and the putative evolution of the hot and dense 
matter created in these collisions into a thermalized, deconfined state of matter called 
a quark gluon plasma. The theoretical challenge is to understand the dynamics of 
the formation of this matter and its properties from QCD while the experimental 
challenge is to detect evidence that such a plasma was indeed formed pj . 

The space-time evolution of the nuclei after the collision and the magnitudes and 
relevance of various proposed signatures of this hot and dense matter depend sensi- 
tively on the initial conditions for the evolution, namely, the parton distributions in 
each of the nuclei prior to the collision. In the standard perturbative QCD approach 
to the problem, observables from the collision may be computed by convolving the 
parton distributions of each nucleus, determined from deep inelastic scattering ex- 
periments, with the elementary parton-parton scattering cross sections. At the high 
energies of the RHIC and LHC colliders, hundreds of mini-jets may be formed in 
the initial collision |3[ B. The final state interactions of these mini-jets are often 
described in multiple scattering Glauber-Gribov models (see Ref. || and references 
therein) or in classical cascade approaches to obtain the space-time evolution (see 
Ref. H and references therein). The possible "quenching" of these mini-jets has 
also been studied and proposed as a signature of the formation of a quark gluon 
plasma |/]]. Recently, initial conditions for the energy density and velocity obtained 
in the mini-jet approach have been used in a simple hydrodynamic model to study 
the late time evolution of matter in high energy nuclear collisions j|, || . 

While the above "probabilistic" approach provides a reasonable description of 
large transverse momentum processes at large x, QCD coherence effects become im- 



portant as we go to small x or alternatively, towards central rapidities [1C]. This 
is because small x partons in one nucleus may "see" more than one parton in the 
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direction of the incoming nucleus resulting in a breakdown of the above described 
convolution of distributions. What is needed therefore to describe the collision of the 
"wee" nuclear partons is a wave picture where coherent multiple scattering is fully 
taken into account. 

At what values of x would this picture become applicable? Strictly speaking, 
coherence or shadowing effects start to become important when wee partons from 
neighbouring nucleons overlap, i.e., when x < l/(2i?jvmjv) ~ 0.1 fm. Here Rn de- 
notes the nucleon radius and tun the nucleon mass. However, several authors have 
shown that with an appropriate choice of the non-perturbative input structure func- 
tions, the standard leading twist perturbative evolution equations still work fairly 
well in the shadowing region below x = 0.1 [11, 12]. To determine at precisely what 
value of x coherence effects become large one needs to compute the two gluon dis- 
tribution function at small x. Following the pioneering work of Gribov, Levin and 



Ryskin [13] and of Mueller and Qiu [14], there has been a considerable body of recent 
work directed towards addressing this issue ]^0|, 15]. 

In this paper, we will describe an ab initio QCD based effective theory approach to 
the theoretical study of nuclear collisions at very high energies. This model of high en- 
ergy nuclear collisions naturally incorporates coherence effects which become impor- 
tant at small x and small transverse momenta while reproducing simultaneously the 
standard mini-jet results at large transverse momenta. It has the further advantage 
of containing a self-consistent space-time picture of the nuclear collison. The model 
is based on an effective action approach to QCD initially developed by McLerran and 
Venugopalan [jl6(|, and later further developed by Ayala, Jalilian-Marian, McLerran 



and Venugopalan [17], and by J. Jalilian-Marian, Kovner, McLerran, Leonidov and 
Weigert g H ||. 

The above mentioned effective action contains one dimensionful parameter, x(y> Q 2 ) 
(This parameter is also often used interchangeably with [i in the text; the distinction 
between the two is discussed in section 3.) Here % is the- total color charge squared 
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per unit area integrated from the rapidity y of interest to the beam rapidity. It is the 
only scale in the problem and we expect therefore that the coupling constant runs as 
a function of this scale. One therefore has weak coupling in the limits where the color 
charge \ is large; either for A S> 1 or s — > oo. It was argued that the classical fields 
corresponding to the saddle point solutions of the effective theory are the non-Abelian 
analogue of the Weizsacker- Williams fields in classical electrodynamics. Exact ana- 
lytical expressions for these fields have been obtained recently [JO], E3J. Further, it 
has been shown explicitly that \ obeys renormalization group equations in y and Q 2 . 
These reduce to the well known BFKL and DGLAP equations respectively |2^] in the 
appropriate limits |l9|, pOp . 

The above model was first applied to the problem of nuclear collisions by Kovner, 
McLerran and Weigert, who formulated the problem as the collision of Weizsacker- 
Williams fields p3|]. The classical fields after the collision then correspond to solutions 
of the Yang-Mills equations in the presence of static, random sources of color charge 
on the light cone. The initial conditions for the dynamical evolution of the small x 
modes of the two nuclei after the collision were formulated and perturbative solutions 
obtained for modes with transverse momenta kf >> «s\/x- After averaging over the 
Gaussian random sources of color charge on the light cone, the energy and number 
distributions of physical gluons were computed. Further, the classical gluon radiation 
from these perturbative modes was studied by these authors and later in greater detail 



by several others [^J, 25, 26]. In the small x limit, it was shown that the classical 
Yang-Mills result agreed with the quantum Bremsstrahlung result of Gunion and 
Bertsch ||. 

While the perturbative approach is very relevant and useful, it is still essential 
to consider the full non-perturbative approach for the following reasons. Firstly, the 
classical gluon radiation computed perturbatively is infrared singular and has to be 
cut-off at some scale. This problem also arises in mini-jet calculations where at high 
energies results are shown to be rather sensitive to the cut-off |J. It was argued in 
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Ref. p3| , 25] that a natural scale where the distributions are cut-off is given by kt ~ 
a S\/X- However, since quantitative differences can be large, it is important to perform 
a full calculation. Secondly, the non-perturbative approach is crucial to study the 
space-time evolution of the nuclei and in particular, the possible thermalization of 
the system and the relevant time scales for thermalization. This in turn has several 
ramifications for computations of various signatures of the quark gluon plasma. For 
instance, if thermalization does occur, then as proposed by Bjorken [27] hydrodynamic 
evolution of the system is reasonable. In that event, our approach would provide 
the initial temperature and velocity profiles necessary for such an evolution || (see 
also ]^8| and references therein). 

We discuss in this paper results from real time simulations of the full, non- 
perturbative evolution of classical non-Abelian Weizsacker-Williams fields. Such a 
simulation is possible since the fields are classical. Similar simulations of the real time 
evolution of classical fields have been performed in the context of sphaleron-mediated 



baryon number violation [30] and chirality violating transitions in hot gauge theories 

ill 



In brief, the idea is as follows [32]. We write down the lattice Hamiltonian which 
describes the evolution of the small x classical gauge fields. It is the Kogut-Susskind 
Hamiltonian in 2+1-dimensions coupled to an adjoint scalar field. For simplicity, 
we restrict our study in this paper to an SU(2) gauge theory. The extension to 
the physical SU(3) case will be considered at a later date. The lattice equations 
of motion for the fields are determined straightforwardly using Hamilton's equations. 
The initial conditions for the evolution are provided by the Weizsacker-Williams fields 
for the nuclei before the collisions. Interestingly, the dependence on the static light 
cone sources does not enter through the Hamiltonian but instead from the initial 
conditions. Also, to reiterate, our results have to be averaged over by the above 
mentioned Gaussian measure for each source. 

A limitation of our approach is that it is classical-quantum fluctuations have been 
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neglected. However if the effective action approach captures the essential physics 
of the small x modes of interest, then in the spirit of the Wilson renormalization 
group, quantum information from the large x modes (above the rapidity of interest) 
is contained in the parameter x{y,Q 2 ) discussed above, which grows rapidly as one 
goes to smaller and smaller x's. This information can be included in the classical 
lattice simulations. Thus as long as we are not at small enough x where the above 
picture of Gaussian random sources breaks down, the residual quantum fluctuations 
about the classical saddle point of the effective theory should be small and the above 
classical picture of nuclear collisions should be valid Q. 

A related approach is that of Mueller, Kovchegov and Wallon p3| , |34| ], where 
they combined Mueller's dipole picture of high energy scattering ^B, 36] with the 



classical Yang-Mills picture [21] to study nucleon-nucleus scattering. In particu- 
lar, Mueller and Kovchegov make the interesting observation in their calculation of 
nucleon-nucleus scattering that, in light cone gauge, the effect of final state interac- 
tions is already contained in the wavefunction of the incoming hadron. This obser- 
vation is very much at the heart of this work because what final state interactions 
there are, are very much determined by the initial conditions given by the small x 
wavefunctions of the nuclei before the collision. For alternative approaches, we refer 



the reader to the work of Makhlin and Surdutovich [37] and that of Balitskii |38| ], 

Our paper is organized as follows. In the following section we discuss the prob- 
lem of initial conditions for nuclear collisions as formulated by Kovner, McLerran and 
Weigert and their perturbative solutions of the Yang-Mills equations and computa- 
tion of gluon production in this approach. We discuss a non-perturbative Hamilto- 
nian approach to the solution of the full Yang-Mills equations. 



The classical picture of nuclear collisions will also be valid at very small x but the Gaussian 
weight with which we average over the classical configurations will change. The functional replacing 
the Gaussian weight is the solution of a non-linear renormalization group equation J2(^] and is yet to 
be determined. 
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In section 3, we formulate the problem of solving the Yang-Mills equations on 
the lattice. Assuming boost invariance and N c = 2, we write down the lattice action 
in 2+1-dimensions. The lattice action is used to construct the lattice analog of 
the continuum initial conditions. This is done by matching the singular pieces of 
the lattice equations of motion on the light cone. Once the initial conditions are 
determined, we write down the lattice Hamiltonian and the equations of motion for 
the evolution of the dynamical fields and their conjugate momenta in the forward 
light cone. 

In section 4, we use lattice perturbation theory to relate the parameters of our 
lattice calculation, such as the parameter proportional to the parton density, fi, the 
lattice size L and the lattice time r, to physical strong interaction scales. We also 
discuss the relation of gauge invariant quantities computed on the lattice to exper- 
imental observables that might be measured in heavy ion collisions at RHIC and 
LHC. 

Numerical results from our simulations are discussed in section 5. These are 
performed for a range of values of g 2 fi = 0.018-0.2, and for lattice sizes from 10 x 10 
to 160 x 160, both measured in units of the lattice spacing. (We will assume boost 
invariance throughout, which simplifies our simulation to a 2-dimensional one.) We 
first compare the field intensities and the time evolution of hard modes in the colli- 
sion with the predictions of lattice perturbation theory and find excellent agreement. 
Significant deviations from lattice perturbation theory are found for the softer modes 
kt ~ as/J-, particularly at larger values of fi. We study the dependence of our re- 
sults on the lattice size. We demonstrate the behaviour of the chromo-electric and 
chromo-magnetic fields as a function of time and that of those components of the 
stress energy tensor that may be related to experimental observables. 

We summarize our results in section 6 and discuss further computations that 
may be performed in this approach. There are two appendices. The first appendix 
discusses our numerical algorithm and procedure. The second is a derivation of the 
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lattice perturbation theory expressions for the field intensity and the kinetic energy at 
r = 0. These expressions were compared to the full lattice results at large transverse 
momenta. 



2 Classical gluon radiation in high energy nuclear colli- 
sions 

In the work of McLerran and Venugopalan [jjql , the classical gluon field at small x 
for a nucleus in the infinite frame is obtained by solving the Yang-Mills equations 
in the presence of a static source of color charge p a (rt,rj) on the light cone. This 
corresponds to the saddle point solution of their small x effective action. Exact 
solutions for the classical field as functions of p a (rt, n) were found by Jalilian- Marian 
et al. [18] and independently by Kovchegov |2l]| . Distribution functions are computed 



by averaging products of the classical fields over a Gaussian measure in p with the 
variance p 2 (rj, Q 2 ). Here p 2 is the color charge squared per unit area per unit rapidity 
resolved at a scale Q 2 by an external probe. It is related to x by the expression 



X (r],Q 2 )= I d7//i 2 (i/,Q 2 ). (1) 

v 



The above picture of gluon fields in a nucleus at small x was extended to describe 
nuclear collisions by Kovner, McLerran and Weigert 23]. 



In this section, we shall discuss their formulation of the problem in the continuum 
and their perturbative computation, to second order in the parameter asp/kt, of 
classical gluon radiation in nuclear collisions. (Readers familiar with the discussion 
in Refs. ^3|-[^(]] can skip sections 2.1 and 2.2 and go directly to 2.3.) In section 
2.3 we then briefly discuss a non-perturbative Hamiltonian approach which suggests 
how all orders in asp/kt can be computed numerically. The implementation of this 
approach on the lattice is described in section 3. 
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2.1 The non-Abelian Weizsacker— Williams approach to high energy 
nuclear collisions 



In nuclear collisions at very high energies, the hard valence parton modes act as highly 
Lorentz contracted, static sources of color charge for the wee parton, Weizsacker- 
Williams modes in the nuclei. The sources are described by the current 



where p\ and p2 correspond to the color charge densities of the hard modes in nucleus 
1 and nucleus 2 respectively. The classical field of two nuclei describing the wee parton 
dynamics is given by the solutions of the Yang-Mills equations in the presence of two 
light cone sources. We have then 



Gluon distributions are simply related to the Fourier transform Af(kt) of the 
solution to the above equation by < Af(kt)Af(kt) > p . The averaging over the classical 
charge distributions is defined by 



The averaging over the color charge distributions is performed independently for each 
nucleus with equal Gaussian weight g 4 p 2 . (Note that this is only true for the case of 
identical nuclei which will be assumed implicitly throughout the rest of this paper.) 

The observant reader will notice that we have omitted the rapidity dependence 
of the the charge distributions in the equations immediately above. We will justify 
this omission in our discussion of the lattice Hamiltonian. We note that the rapidity 
dependence of the charge distribution is also absent in Ref. |2^] (see the discussion 



r> a {r t ) = 5 v +pl{r t )8{x-) + V-(j&n)5{x+) , 



D^F^ = J' 



(3) 




below Eq. 0). 



S 



Before the nuclei collide (t < 0) , a solution of the equations of motion is 
A ± = 0, 

A i = 6(x-)6(-x + )a[(r t )+8(x + )e(-x-)a 2 (rt), (5) 

where a q (r t ) (q = 1,2 denote the labels of the nuclei) are pure gauge fields defined 
through the gauge transformation parameters A q (rj,rt) p4| ] 

a q (r t ) = -[Pe J± ^ ' V Pe J± Vi q(V *M . (6) 

Here 77 = 7? pro j — log(x~/Xp ro j) is the rapidity of the nucleus moving along the positive 
light cone with the gluon field a\ and 77 = — r/ pro j +log(Xp ro j/2; + ) is the rapidity of the 
nucleus moving along the negative light cone with the gluon field a\- The A q (rj,rt) 
parameters are, in turn, determined by the color charge distributions: 

A ± A q = Pq ■ q = l,2 (7) 

being the Laplacian in the perpendicular plane. 

It is expected that at central rapidities (or x <C 1) the source density varies 
slowly as a function of rapidity and a 1 = a l (rt). The above expression suggests that 
for t < the solution is simply the sum of two disconnected pure gauges. Just as 
in the Weizsacker-Williams limit in QED, the transverse components of the electric 
field are highly singular. 

For t > the solution is no longer pure gauge. Working in the Schwinger gauge 

x + A~ +x~A + = 0, (8) 
or A T = 0, the authors of Ref. |^3[ found that with the ansatz 

A^ = ±x ± a(r, r t ) , 

A* = a l ± (r,r t ), (9) 
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where r = V2x + x , Eq. could be written in the simpler form 



rd T T 3 d T a + [Di, 



D\a 



1 



[Di, d T a\} + ir[a, d T a] 



1 



d T Td T a l > 



ir 2 \a, 



D\a}- [D j ,F ji 



0. 
0. 
0. 



(10) 



Note that the above equations of motion are independent of rj-the gauge fields in the 
forward light cone are therefore only functions of r and rt and are explicitly boost 
invariant. We will use this fact later in our discussion of the Hamiltonian approach. 

The initial conditions for the fields a(r, rt) and a\ at r = are obtained by 
matching the equations of motion (Eq. |3|) at the point = and along the bound- 
aries x + = 0, x~ > and x~ = 0, x + > 0. Because the sources are highly singular 
functions along their respective light cones, so too in general will be the equations 
of motion. Remarkably, however, there exists a set of non-singular initial conditions 
that ensure the smooth evolution of the classical fields in the forward light cone. 
These are obtained by matching the singular terms in the equations of motion before 
and after the collision at r = 0. In terms of the fields of each of the nuclei before the 
collision (t < 0), the fields after the collision (t > 0) are 



a±\ T =o 
a\ T =o 



(11) 



Gyulassy and McLerran have shown [gj] that even when the fields a\ 2 before the 
collision are smeared out in rapidity to properly account for singular contact terms in 
the equations of motion the above boundary conditions remain unchanged. Further, 
since the equations are very singular at r = 0, the only condition on the derivatives 
of the fields that would lead to regular solutions are d T a\ T= o, d T a\\ T= Q = 0. 
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2.2 Review of perturbative solution of the Yang— Mills equations 

In Ref. [^jj , perturbative solutions (for small p) were found to order p 2 by expanding 
the initial conditions and the fields in powers p (or equivalently, in powers of asp/kt) 
as 

oo oo 

a = £ a (n) ; a ± = 2 a ^(«) ' ( 12 ) 

n=0 n=0 

where the subscript n denotes the nth order in p. Since it is relevant to the discussion 
in the following sections (particularly appendix B) , we outline their solution below 
with a commentary but refer the reader to their paper for the details. 

To lowest order in p, the Yang-Mills equations of motion are linear in am and 
a ±(i)' som tion is the trivial Abelian solution 

«(i) = ; = - + ct> 2 )(x t ) , (13) 

where (f) q = ^rP q - Clearly aj_(i) above is a pure gauge. To second order in p, the 
equations of motion are nearly homogeneous except for a non-linear term proportional 
to a J_(i)] m ^ ne equation for CKj_/ 2 \. This piece can be removed from the 

equations of motion to this order by making use of the residual, r-independent gauge 
transformation = V(xt)[e fl — jd^] V'(xt) and fixing V such that e l satisfies the 
two dimensional Coulomb gauge condition d l e l \ T= o = 0. As described in appendix 
B, the lattice Coulomb gauge condition is fixed in an analogous way to eliminate the 
residual gauge freedom on the lattice. 

Writing e l = e^d^Xi where e lJ is the two dimensional anti-symmetric Levi-Civita 
tensor, the Yang-Mills equations in the forward light cone (Eq. |l0|) may be written 
as 

-L<9 T T 3 d r e ( 2) - Vj.e(2) = , 

-d T Td TX{ 2) - Vix ( 2) = 0. (14) 
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with the initial conditions 



e (2)|r=0 
X(2)|r=0 



-\d^\d^ 2 }, 



(15) 



The residual Coulomb gauge fixing at this second order therefore shifts the non- 
linearities from the equations of motion to the initial conditions. 



The solutions of Eqs. 14 can be written in terms of Bessel functions as 

d 2 k t d 2 y t 



e (2){T, X t ) 
X(2)(l~, X t ) 



(2vr) 2 
(2tt)2 



^-^(ift) — Ji(wr), 



•jJT 



(16) 



where 



h 3 (yt) 
h\{yt) 



1 



V 



(17) 



Here ui = \kt\- Note that in the solutions to the gauge fields above, the spatial and 
temporal distributions factorize. The amplitudes of the fields in this perturbative 
approach are therefore completely determined at r = 0. In section 5, we will see that 
in weak coupling the large transverse momentum modes show the above Bessel be- 
havior. (In lattice perturbation theory, the form is the same as above. The difference 
is that uj obeys a lattice dispersion relation). 

At late times ujt >> 1, the well known asymptotic expressions for the Bessel 
functions can be used to write the solutions in Eq. 16 as 

d 2 k t 1 



J(r,x t ) = 
Here k % = e % ^k\juj and 
axikt) 



(2tt) 2 

d 2 k t ij_ 
(2vr) 2 " 

7T UJ 



ai(k t 



a 2 (k 



1 

1 



ikfXt—ioiT 



+ h.c 



-1/2 



ikfXt—iuiT 



+ h.c 



; a 2 (k t ) = —=iujhi(k t ) , 



(18) 



(19) 
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where hi and /13 are now the Fourier transforms of Eq. 17. 

The energy distribution in a transverse box of size R and longitudinal extent 
dz can be computed by summing over the energy of the modes in the box with the 
occupation number of the modes given by the functions a,i(kt) above. We have then 
(for cut » 1) 

The multiplicity distribution of classical gluons is defined as dE/dyd 2 kt/u>. After per- 



forming the averaging over the Gaussian sources, the number distribution of classical 
gluons is 



dyd?k t (2vr) 4 kf 

where L(kt,\) is an infrared divergent function at the scale A. It will be discussed 
further below. We first note that this result agrees with the quantum bremsstrahlung 



formula of Gunion and Bertsch [29] and with several later works [24, [26|]. It was 
also shown by Gyulassy and McLerran that when the sources are smeared in rapidity, 
the expression that results is identical to the one above except /i 4 — ► X + {y)x iy) 
where the ± superscripts refer to the nucleus on the positive or negative light cone 
respectively. 

The origin of the infrared divergent function L(kt,X) above is from long range 
color correlations which are cut-off either by a nuclear form factor (as in Refs. [29|, 25[|), 



by dynamical screening effects |39| , |40| or in the classical Yang-Mills case of Ref . [P3| 1 , 
non-linearities that become large at the scale kt ~ as/J. In the classical case then, 

L(^,A)=log(^/A 2 ), (22) 

where A = ag/x. A similar logarithmic behaviour at small transverse momenta can 



be deduced from the dipole form factors used in Refs. [29, 39]. A power law (~ 1/A;, 



infrared behaviour is predicted in Ref. . The formalism used in all these derivations 
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breaks down at small momenta and one cannot distinguish between the different 
parametrizations of the nuclear form factors. However, at sufficiently high energies, 
the behaviour of L(kt,X) in the infrared is given by higher order (in asf^/kt) non- 
linear terms in the classical effective theory. One of the goals of our work is to address 
precisely this question: how do non-perturbative effects in the classical effective 
theory change the gluon distributions at small transverse momenta? 

2.3 The Hamiltonian approach 

While the Yang-Mills equations can be solved perturbatively, in the limit ctsfJ- -C kt, 
it is unlikely that a simple analytical solution exists for Eq. || in the non-perturbative 
regime where kt < /U. The classical solutions have to be determined numerically for 
t > 0. The straightforward procedure would be to discretize Eq. ^. It will be more 
convenient for our purposes though to construct the lattice Hamiltonian and obtain 
the lattice equations of motion from Hamilton's equations. This will be done in the 
next section. Before we do that, we will discuss here the form of the continuum 
Hamiltonian and comment on our assumption of boost invariance. 
We start from the QCD action (without dynamical quarks) 

Sqcd = J d^x^g-^-g^g^F^ - f , (23) 

where g = det(gw). In the forward light cone (t > 0) it is convenient to work with the 
r,r],ft co-ordinates where r = V2x + x~ is the proper time, rj = \ log(x + /x _ ) is the 
space-time rapidity and ft = (x,y) are the two transverse Euclidean co-ordinates. 
In these co-ordinates, the metric is diagonal with g TT = —g xx = —g vy = 1 and 
g vv = -1/r 2 . 

After a little algebra, the Hamiltonian can be written as [JO]] 
H = rJ dT^ t {ij^' a + ^P r ^ -(24) 
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Here we have adopted the gauge condition of Eq. ||, which is equivalent to requiring 
A T = 0. Also, p 71 = ^pdrAjj and p r = rd T A r are the conjugate momenta. 

Consider the field strength in the above Hamiltonian. If we assume approx- 
imate boost invariance, or 

A r (T,rj,r t ) » Ar(r,fi)\ A v (T,rj,7%) ss $(<r,r t ), (25) 

we obtain 

F%r = -D r $ a , (26) 

where D r = d r — igA r is the covariant derivative. Further, if we express j v ' r in terms 
of the defined in Eq. ^ we obtain the result that j^ ,r = for r > 0. 

Performing the integration over the space-time rapidity, we can re-write the 



Hamiltonian in Eq. 24 as 

H = j dr^^E^E* + T -F? 5 F% + i- (Z^r (D r <S>f + Ip^p"} . (27) 



We have thus succeeded in expressing the Hamiltonian in Eq. 24 as the Yang-Mills 



Hamiltonian in 2+1-dimensions coupled to an adjoint scalar. The discrete version of 



the above Hamiltonian is well known and is the Kogut-Susskind Hamiltonian [42] in 
2+1-dimensions coupled to an adjoint scalar field. The lattice Hamiltonian will be 
discussed further in the next section. 

We now comment on a key assumption in the above derivation, namely, the boost 



invariance of the fields. This invariance results in Eq. 27 thereby allowing us to restrict 



ourselves to a transverse lattice alone. To clarify the issue we are compelled to make 



a few historical remarks. As we mentioned earlier, the authors of Ref. [23] found a 
solution which was explicitly boost invariant. However, this result was a consequence 
of the original assumption of McLerran and Venugopalan that the color charge density 



factorizes, p a (rt, rj) — > p a (rt)5(x ). It was noticed in Ref. 43] that this factorized form 



for the charge density results in infrared singular correlation functions which diverge 
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as the square of the lattice size. This problem was subsequently resolved in Ref. [18 
where the authors realized that a rapidity dependent charge density p a (rt,r/) would 
give infrared safe solutions. This might be interpreted as implying that the boost 
invariance assumption of Ref. 23] should be given up as well. 

Fortunately, this is not necessary. In principle, the rapidity dependence of the 
color charge density can be arbitrarily weak since that is sufficient to obtain infrared 
safe correlation functions. In Ref. [44|, an explicit model was constructed for the 
color charge distribution in the fragmentation region. It was shown there that for 
Tj < f?proj the color charge distribution had a very weak dependence on rj. Further, 
as mentioned earlier, it was shown by Gyulassy and McLerran |24|] that the initial 
conditions in Eq. |ll] are unaffected by the smearing in rapidity. 

In general, at the energies of interest, particle distributions are unlikely to be 
boost invariant. Eskola, Kajantie and Ruuskanen M have shown in the mini-jet 
picture that the final distributions are more like broad Gaussians. This would be 
true in our case as well since fi 2 is in general a function of the rapidity and this 
dependence may be strong in the central region. If we wish to describe particle 
distributions for a range of rapidities as opposed to our current restriction to 1 unit 
in rapidity, we will have to give up our assumption of boost invariance. This can be 
easily done, but it will be numerically more time consuming as well. 



3 Real-time lattice description of nuclear collisions 

In the previous section we discussed the continuum formulation of nuclear collisions in 
terms of collisions of non-Abelian Weizsacker-Williams fields. The initial conditions 
for the evolution after the collision were obtained by requiring that the singular con- 
tributions to the equations of motion on the light cone vanish. We next discussed the 
perturbative solutions of the equations of motion for our choice of initial conditions 
and obtained a result for the number distribution of radiated gluons. 
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In this section, we will formulate the problem on the lattice. We begin by writing 
down the lattice action and equations of motion. From the lattice action, we can 
indentify the singular terms in the lattice equations of motion. Matching these on 
the light cone uniquely gives us the initial conditions on the lattice. We end with a 
discussion of Hamilton's equations of motion for the evolution of dynamical fields on 
the lattice with the initial conditions specified at r = 0. 

3.1 Lattice action and equations of motion 

In this sub-section, we will write down the lattice action in 2+1-dimensions and 
the rules to derive the lattice equations of motion from this action. In the following 
sub-section, we will identify the singular terms in the lattice action and use these to 
derive the correct initial conditions for dynamical evolution of fields on the lattice for 
r > 0. 

The action is defined in the 4-space discretized in the transverse directions, while 
z and t are continuous. The appropriate expression is derived starting from the 
Minkowski Wilson action in the discretized 4-space and taking the naive continuum 
limit in the longitudinal directions. The Minkowski Wilson action for the SU(N C ) 
gauge group in fundamental representation reads 

s = «- 2 E(i-^™*)+E(i-^™.x) 

where zt, z _L, t _L and _L are, in obvious notation, plaquettes lying in various 2-planes 
of the 4-space and 3? denotes the real component. A plaquette is defined as 

Ujim = U jt iU j+hm U j+m l U j m . 

where j is a site index and l,m are direction indices. We now take the formal 
continuum limit in the longitudinal directions by writing longitudinal links as U = 
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exp(iaA), letting a — > and assuming that all the fields are smooth functions of 
longitudinal coordinates. The powers of a in front of various terms in S have been 
chosen with this formal limit in mind. Replacing a 2 ^2 zt with J dzdt, we then have 
for the action 



S = J dzdtY, 



' TrF^ t + -^-5R Tr(M t j_ - M zl _) - ( 1 - -^TrCTj 



2N C Zl N c v - V N c 



(28) 



where 

1 



l -d 2 t ul n + i(A ttj+n d t ul n - d t ul n A td ) + A t , j+n ul n A td 

(29) " 

and similarly for M z j n . 

The equation of motion for a field is obtained by varying S with respect to that 
field. For the longitudinal fields At )Z the variation has the usual meaning of a partial 
derivative. For transverse link matrices U± the variation amounts to a covariant 
derivative 

D~<S{U ± ) = d r S {exp(ira^)U ± ) | r=0 , (30) 

where <t 7 , 1 < 7 < N 2 — 1 form a basis of SU(N C ) algebra.. In particular, 

D~ f U± = ia^U ± ; D^u[ = -iU[a y , (31) 

and derivatives of more complicated functions are obtained, combining these two rules 
with the usual rules of differentiation. 

3.2 Initial conditions on the lattice 



We now derive the lattice analogue of the continuum initial conditions in Eq. 11. We 



start from the lattice action in Eq. 28 , obtain the lattice equations of motion in the 
four light cone regions and determine non-singular initial conditions by matching at 
r = the coefficients of the most singular terms in the equations of motion. 
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On the lattice, the initial conditions are the constraints on the longitudinal gauge 
potential and the transverse link matrices U± at r = 0. The longitudinal gauge 
potentials are zero outside the light cone and satisfy the Schwinger gauge condition 
(Eq. ||) inside the light cone x± > 0. Thus they can be written as in the continuum 
case (see Eq. ||) as 

A ± = ±x ± 9{x + )0{x~)a(T, x t ) . (32) 

From our discussion in the last section, it is clear that the transverse link matrices 
are, for each nucleus, pure gauges before the collision. This fact is reflected by writing 

U ± = 6{-x + )e(- X -)I + 9{x + )9{x-)U(t) + e(-x + )6{x-)U {1) + 9(x + )9{-x~)U (2) , 

(33) 

where C/Vm 2 ) are pure gauge. 

The pure gauges are defined on the lattice as follows. To each lattice site j we 
assign two SU(iV c ) matrices Vij and V-zj- Each of these two defines a pure gauge 
lattice gauge configuration with the link variables 

U$ = V q>j vl j+n , (34) 

where q = 1,2 labels the two nuclei. As in the continuum, the gauge transforma- 
tion matrices V q j are determined by the color charge distribution p q j of the nuclei, 
normally distributed with the standard deviation g^fi^ (compare with the continuum 
distribution in Eq. ||): 

P [Pd] « ex P |"^|E Pg,j j • ( 35 ) 

Parametrizing V q j as exp(iA|) with Hermitean traceless A|, we then obtain A^ by 
solving the lattice Poisson equation 

A L AJ = £ ( A] +n + A]_ n - 2 A]) = Pq>j . (36) 
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It is easy to verify that the correct continuum solution (Eqs. |5]and|9|) for the transverse 
fields A± is recovered by taking the formal continuum limit of Eq. 53l Using the 



general representation of the gauge fields in Eqs. 32 and 33, we shall now derive the 
initial conditions for them at r = 0. 

Initial conditions for x ± > 0: U± 

The equation of motion for U±, contains, upon substitution of U± from (^) and j4 ± 
from (|32"D, two types of terms. 



1. Terms regular at = 0. All the terms containing less than two longitudinal 
derivatives belong to this category. The regularity of terms containing no longi- 
tudinal derivatives is clear. Terms in the action containing a single longitudinal 
derivative all involve the combination A~d- + A + d+. Given the functional form 
(|32|) of j4 ± , these terms only give rise to regular contributions to the equation 



of motion. The double-derivative contributions WTrU±d+d-U± in the action 
give rise to terms in the equation of motion behaving as 

-x ± 5(x ± )9(x T )d T U , 
r 

which is regular if d T U vanishes rapidly enough as r — ► 0. The terms regular 
at x± = provide no relation between U± and U^'^ 2 \ 

2. The singular terms containing the product 5(x + )5(x~). These originate in 
the double-derivative contributions in the action, when both 

derivative operators act on the step functions. Since the coefficient in front of 
5(x + )5(x~) must vanish in order to satisfy the equation of motion, a matching 
relation between U and CT^ 1 )^ 2 ) is obtained. 

We now compute the singular part of the U± equation of motion. We note that 

U±(d? - d 2 z )u[ = 2U±(d+d-U[) = 2(d + d-U ± )u[ 
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modulo a total derivative and utilize the identity 9(x)5(x) = \5{x) with the rules 
of Eq. |3l[ This gives, upon substitution of (Eq. |33| ) and retaining only the terms 
proportional to 5(x + )S(x~), 



Tr(j 7 



(rf. 1 ) +U®)(I + rf)-h.c] = 0. (37) 



We therefore have obtained the result that (U^ + U^)(I + U^) should have no anti- 
Hermitean traceless part. Note that this condition has the correct formal continuum 
limit: writing [/( 1 )'( 2 ) as exp^a^a^) and U as exp(ia±a±), we have, for small a±, 

a± = a\ + 02, 

as required. 

The above condition in Eq. ( |37|) can easily be resolved in the SU(2) case, because 
the sum of any two 2x2 unitary matrices is a unitary matrix times a real number. 
Also, the anti-Hermitean part of an SU(2) matrix is traceless, while the Hermitean 
part is proportional to I. Therefore in this case Eq. [37] is equivalent to 

(UW + uW)(i + tf) = ci, 

where C is a real number. Resolving this equation for W and requiring that UW = I, 
we obtain, after a simple calculation, that C = Tt(U^ + U^) and hence 

u = T ^ uil) + U(2) ) _ j = {u m + + rfnyi . (38) 

1/(1)1 + ^( 2 )T 

For N c > 2 the solution can be obtained by solving the — 1 equations (^) for the 
— 1 real numbers parametrizing U. The solution does not necessarily have the 
compact form ( |38|) . For simplicity, we leave the N c > 2 case outside the scope of this 
paper. 

Initial conditions for x ± > 0: A ± 

We have seen that the matching conditions for the fields before and after the collision 
stem from the requirement that the singular terms in the equations of motion cancel 
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out. We now impose this requirement on the equation of motion for A_. As before, 
the singularities arise only because longitudinal derivatives of step functions are taken. 

Consider first the i* 1 ? term in the action. Since in this case we need two deriva- 
tives for a genuine singularity, we are only interested in the Abelian part of F?_, 
whose variation with respect to A +, ' y gives 

1 



N c 



Tia 1 d + (d-A + - d+A„), 



whose most singular part is (using Eq. 32 and x5'{x) 

a^6(x~)5(x + ) . 



-S(x)) 



We now vary the ±, _L terms (Eq. |29| ) in the action (Eq. |28|) with respect to A - ' 7 and 
select the contributions containing derivatives. The result is 

]TTra 7 [(d + U]_ njn )U^ n>n - U^d+Uj^) - I ].»{<)■ I j.„ ! + (d+U^U^ 



2N< 



the singular part of which is 



+ {U^~ C/ (2) ),-„,n-(C/ (2)t -C/ (2) kn • 
Requiring that the coefficient in front of 5{x + )9{x~) vanish, we obtain 

^ J2 Tra, {(U^ + UO-W - h.c.) i|B - (U^ + rfV® - h.a)i 



(39) 



This result can be cast in a more symmetric form, using the initial condition for U±. 
Note that 

C/(2) + u{i)u\ = i (C/ (i) + [/(2) )(/ + ^t) + _ £/(2))([/t _ j) . 

Since the first term on the right-hand side has no anti-Hermitean traceless part, it 
follows that 

i 



a 



7 



4/V, 



^TrcT 7 ([(C/W - uW)(tf -I)- h.c.]i,„ - [(f/ f - /)(£/« - C/ (2) ) - h.c/ 
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It is easily seen that the above equation has the correct formal continuum limit. 
Writing again f/Wil 2 ) as exp^a^a^) and U as exp(ia±a±), we have, for small a±, 

Hence, in the limit of smooth fields, 

a = i^2[ai,a 2 ] n , 

n 

as required. 



Finally, we require that the last line of Eq. 39 be cancelled by the term in the 
action coupling A + to the external source. This term then must have the form 

JL / dzdt$>(z+)TrA+£ \{U^ - U^ n>n - (U^ - UW) j>n , 
whose formal continuum limit is easily seen to be the correct one: 

~w I dAx6 ( x+ ) TtA+v ■ A ± ■ 

A completely analogous term should be included for A~ . 
3.3 Hamiltonian formulation on the lattice 

In the previous sub-sections, we wrote down the lattice action and lattice equations 
of motion for all regions of the light cone. Matching the singular terms on the light 
cone in section 3.2, we obtained the initial conditions for evolution in the forward 
light cone. 

In section 2.3, we derived the continuum Hamiltonian for the forward light cone 
(t > 0) to be the Yang-Mills Hamiltonian in 2+1-dimensions in the gauge A T = 0. 



The lattice Hamiltonian is obtained by performing a Legendre transform of Eq. 28 



following the standard Kogut-Susskind procedure [42]. The analog of the Kogut- 
Susskind Hamiltonian here is 

l=(j,h) D 
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TV 

Figure 1: The lattice size dependence of the scalar kinetic energy density, expressed 
in units of /i 4 for fj, = 0.0177 (pluses) and fi = 0.035 (diamonds). The solid line is the 
LPTh prediction. The error bars are smaller than the plotting symbols. 



4r ^ 

3,n 



(41) 



where E\ are generators of right covariant derivatives on the group and Uj fi is a 
component of the usual SU(2) matrices corresponding to a link from the site j in the 
direction h. The first two terms correspond to the contributions to the Hamiltonian 
from the chromoelectric and chromomagnetic field strengths respectively. In the last 
equation, $ = $ a <j a is the adjoint scalar field with its conjugate momentum p = p a a a . 
Taking the continuum limit of the above Hamiltonian, one recovers the continuum 



Hamiltonian in Eq. 27. 



Lattice equations of motion follow directly from Hl of Eq. [0]. For any dynamical 
variable v with no explicit time dependence v = {Hl,v}, where v is the derivative 
with respect to r, and {} denote Poisson brackets. We take Ei, Ui, pj, and <&j as 
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independent dynamical variables, whose only nonvanishing Poisson brackets are 

{Pi, *$} = SijSab; {Ef, U m } = -i5 lm U l( T a ; {Ef, E b m } = 25 lm e abc Ef 

(no summing of repeated indices). The equations of motion are consistent with a set 
of local constraints (Gauss' laws). These are 
1 



E 



ESfiTr [o a U jfi o>U* A ) - El 



2e abcP )®1 = 0. 



(42) 



Explicitly, using the Poisson brackets above, the equations of motion for the four 
dynamical variables are as follows: 



E, 

in 



-U m Ei 



{ri;(i-i™b),«} 



1 

T 



2$, 



(43) 

(44) 
(45) 
(46) 



Above, $ = U j)fi W] fi . 

The results of section 3 can be summarized as follows. The four independent 
dynamical variables are Ei, U±, pj and <&j. Their evolution in r after the nuclear 
collision is determined by Hamilton's equations above and their values at the initial 
time r = are specified by the initial conditions derived explicitly in section 3.2. For 
easy reference, these are written compactly as follows: 



U\ T =0 
Pj\r=0 



{U X + U 2 ){U\ + Ut)- 1 ; ^|r=0 = . 
2a; $, =0, 



(47) 



where U and a are given by Eq. 38 and Eq. respectively. Note that the second 



set of conditions for $j and pj follows respectively from Eq. |25| and the definition of 



Pj-see the discussion after Eq. 24. 
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4 Interpreting lattice results for continuum physics 

Before we turn to discussing in detail the numerical results from our lattice simula- 
tions, we should consider first the ramifications for the problem of nuclear collisions 
at very high energies. Specifically, we wish to know how one interprets in physical 
terms the results of our numerical simulations. 

In section 2, we introduced a scale fi 2 , the color charge squared per unit area, and 
argued that new non-perturbative physics arises at momenta of order kt ~ ctsU- In 
the original work of McLerran and Venugopalan, only the valence quarks were taken 
to be sources of color charge which gave fi 2 ~ A 1 / 3 fm~ 2 . Hence [i 3> Aqcd only 
for nuclei much larger than physical nuclei. However, if semi-hard gluons at x values 
greater than those of interest here f] are also included as sources of color charge as 
they should be, then fi 2 is significantly larger. Then fi 2 is defined as ||24f 



where q, g stand for the nucleon quark and gluon structure functions at the resolution 
scale Q of the physical process of interest. Also, above xq = Qj \fs. Using the HERA 
structure function data, Gyulassy and McLerran estimated that /U < 1 GeV for LHC 
energies and \x < 0.5 GeV at RHIC. Thus the regime where the classical Yang- 
Mills picture can be applied is rather limited. However, a window does exist and 
depending on what higher order calculations will tell us, this window may be larger 
or smaller than the naive classical extimates. Moreover, the insights gained from 
this self-consistent spacetime approach may still be very valuable in modelling heavy 
ion collisions which are not at the asymptotic energies where the approach becomes 
quantitative. 

This work is not a quantitative study to be compared to experiment. Rather, 
2 Note that in this case, these correspond to those values of x corresponding to rapidities greater 
than the central rapidity in the nuclear collision. 




(48) 
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our goal is a qualitative understanding of non-perturbative phenomena in the central 
region of a heavy ion collision. Towards this end, we need to 

• relate lattice parameters to physical scales, 

• define and compute quantities on the lattice that may be eventually compared 
to experiments with as little ambiguity as possible. 

Let us first relate lattice units and physical units. Given A of a nucleus, we should 
match valence parton densities. To this end, we require the equality of cross-sectional 
areas. The relation between the linear size L = Na of the lattice and the transverse 
radius R of the nucleus is then L 2 = ttR 2 . We ignore for the moment the difference 
in the lattice and continuum boundary conditions. For an A = 200 nucleus R = 6.55 
fm translates into L = 11.6 fm. We can then express the lattice spacing a = L/N in 
units of fermi. 




0.5 1 1.5 2 2.5 

T 

Figure 2: Normalized field intensity of a hard (k t = 2.16GeV) mode vs proper time 
r in units of fm (diamonds). Solid line is the LPTh prediction. 

Next, we determine the relation between the continuum fi and its lattice coun- 
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terpart. To do so, we use the Poisson equation (^) for the gauge transformation 
parameter A: A^A = pi, and compare it to the continuum equation ([?]): Aj_A = p c . 
For a — ► the lattice Laplacian approaches a 2 A, hence we must have pi = Pc® 2 - 
Also, as a — > 0, the lattice and continuum distributions of the color charge should be- 
come identical. Keeping in mind that a 2 J2 = J dx^ for small a, we obtain pi = ap c . 



In the effective theory on the lattice the coupling is g 2 p^L [43|. In this qualitative 
work, for simplicity, we choose g = 1. Then pi = 0.1, for instance, for an A = 200 
nucleus corresponds to 

p c = — 0.1 — ~ 0.28 GeV for N=160 . (49) 

For a physical scale of say p = 0.6 GeV, and the number of sites N = 160, we should 
pick pi ~ 0.21. To compare lattice momenta to momenta in physical units, note that 
h = 2vrn/L, where —{N - l)/2 < n < (N - l)/2. Then, for L = 11.6 fm, k t ~ 1 
GeV for n ~ 9-10. For p c = 1 GeV or equivalently, /i£, = 0.36, one may expect 
perturbation theory to be reliable for n > 10. One can in principle estimate the value 
of n where strong coupling effects become significant. 

There still remains the question of how to relate lattice time ri a tt to the time in 
physical units, T p h ys . This can be obtained by comparing the continuum and lattice 
Hamiltonians at small a, [ 30 1 



Tphys "■''"lattice • 

Recall that a = L/N fm. For L = 11.6 and N = 160, a = 0.07 fm. 

Now that we have straightened out the relation of lattice units to physical units, 
we must consider what can be measured on the lattice. Interesting quantities would be 
those for which our simulations would predict non- perturbative corrections to "em- 
pirical" quantities which can be computed in perturbative QCD at large momenta. 
A quantity one may consider is the cross section for gluon mini-jet production. At 



large transverse momenta, Gyulassy and McLerran [24| have shown that the clas- 
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sical Yang-Mills formula in Eq. 21 is at small x (approximately) the same as the 



perturbative QCD prediction 13] for the process AA — > g 



do a s N c f 2 f(x 1 ,q 2 )f(x 2 ,(k t -q t ) 2 ) 



dyd 2 k t 7T 2 k t 2 J Ht q 2 (k t -q t ) 2 

where 

f^-^w 1 - < 51 > 

and x\ ~ X2 = kt/y/s. The two formulae are equivalent if we divide the above formula 
by 7ri? 2 , approximate the integral above by factoring out / above at the scale k 2 and 
taking the normalization factor Kn ~ 5. At large transverse momenta therefore, we 
can relate the field intensity measured on the lattice |^4fcJ 2 for the appropriate values 
of n to the cross section for AA — ► g. The lattice perturbation theory expression 



(which matches the prediction of the full theory at large momenta) is given by Eq. 39 
of appendix B. 

However, at smaller transverse momenta, there is no simple relation between the 
field intensity and the classical gluon distribution function (and thereby the cross 
section by the above arguments). This is because the dispersion relation for soft 
gluons is not that for free gluons-the interpretation of the field intensity as the gluon 
distribution is then clearly not right. We have to look for more general gauge invariant 
quantities which, conversely, in the limit of large kt, will give us the AA — * g cross 
section. 

Some of the quantities we can compute on the lattice in strong coupling and in 
principle compare to experiment are quantities related to various components of the 
energy-momentum tensor. In an interesting recent paper, Testa has shown |53|] that 
a class of semi-inclusive observables (such as, for example, the energy flux through 
a surface corresponding to the resolution of a detector) can be related to various 
components of the energy-momentum tensor. Our eventual goal is to compute in the 
classical effective theory, non-perturbative corrections to the energy-energy correla- 



tors measured in jet physics [54] 
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Figure 3: Field intensity over fi 4 as a function of kt for \i = 200MeV (squares), 
= lOOMeV (pluses), and = 50MeV (diamonds). Solid line is the LPTh prediction. 
The field intensity is in arbitrary units and kt is in GeV. 

In this work, we measure a variety of observables. On the one hand, extensive 
quantities like the energy, with well defined local gauge-invariant densities, provide 
useful information on the system and are easily measurable. On the other hand, in 
an experiment one usually deals with gluon fields with definite momenta. Such fields 
are not easy to define outside perturbation theory. In order to study momentum 
distributions of gluons, we introduce two types of quantities. 

First, we define the gauge field on every link of the lattice as the imaginary part 
of the link matrix. Note that we work in the A T = gauge and fix the Coulomb 
gauge in the transverse plane for r = 0. Hence there is no gauge freedom left in our 
definition of the gauge field. The same applies for the scalar field. In the following, 
we will study the momentum dependence of these gauge-fixed objects. 

Secondly, we consider a gauge-invariant estimate of the gluon energy in terms of 
equal time energy-energy correlators. We make use of the fact that each of the kinetic 
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and potential terms in ( |4l| ) are squared quantities such as E 2 and B 2 . Consider, for 
instance, the connected correlator in the Abelian case 

C m (x t ) = (B 2 (x t )B 2 (0)) p - (B 2 (x t )) p (B 2 (0)) p . (52) 

If the distribution of B were normal, P[B] oc exp(h[B]), with a translation-invariant 
bilinear functional h, C m (xt) would factorize: 

C m (x t ) = 2{B(x t )B(0)) 2 p > 0. (53) 

Following the analogy with this special Abelian case, we define the magnetic energy 
of the momentum mode kt in the non-Abelian theory as 

E m (k t ) ^^E eikt ' x V^)- ( 54 ) 

The analog of ( |54| ) can be defined for each of the six potential and kinetic terms 
(counting the chromo-electric and the scalar potential terms separately for each lattice 
direction) in (|4l|). The total energy distribution in momentum space is then defined 
as the sum of these six contributions. 

Before we end this section, we should comment on the continuum limit in this 
approach. The theory contains dimensionful quantities, [i and L, which may be 
related to physical observables. One approaches the continuum limit by keeping 
fiL fixed and taking fi in lattice units to zero. Comparing dimensionless ratios of 
physical observables to their lattice counterparts would then allow one to extract the 
continuum value of \xL. In this theory, the field amplitudes squared fall off as 1/k 4 
as opposed to a thermal field theory where the field amplitudes squared fall off as 
1/k 2 . It is therefore more likely that a wider range of physical observables will have 
a well defined continuum limit than in the thermal case. Whether this is indeed the 
case requires a careful study of physical (gauge invariant) observables by taking the 
continuum limit of their lattice counterparts in the manner prescribed above. We 
plan to continue this study in a future work. 
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Finally, we mention here for future reference a proper-time lattice computation 
of an interesting gauge- invariant quantity that can be directly related to experiment, 
namely, the Poynting vector. The Poynting vector describes the energy flux out of 
the nuclear surface. It, and other components of the energy-momentum tensor, can 
be defined on the lattice Computing the Poynting vector however requires a 

more careful choice of boundary conditions than the periodic boundary conditions 
described in this work. This too will be left to a future study. 

5 Real time lattice computation of gluon production: re- 
sults 

We now turn to numerical results from our simulations. These were performed for 
a variety of lattice sizes L= 20-160 and the color charge density /i = 0.015-0.2 in 
units of the lattice spacing a. To convert lattice results to physical units, we take 
L 2 ~ ttR 2 for ^4=200 nuclei. This then also determines ji for a fixed lattice size. The 
relation of lattice time to continuum time is given by the relation tq = a ■ tl [|30|]. 
In Fig. H|, we plot the Gaussian averaged initial kinetic energy (Ek) on the lattice as 
a function of the lattice size L and compare it with the lattice perturbation theory 
expression (Eq. |73| in appendix B) for different values fJ-L = 0.0177, 0.035 of the color 
charge density. For small values of L, there is very good agreement between the two 
but for the largest value L = 160, they begin to deviate. Since the strong coupling 
parameter on the lattice is oc g 2 ^L, for g^^^L > 1 we can expect to see deviations 
from lattice perturbation theory. 

In Fig. H], we plot the ratio of the field intensity of a particular mode of the 
transverse gauge field as a function of proper time t normalized to its value at r = 0. 
The diamonds are results from a lattice simulation with L = 160 and hl = 0.0177 and 
the mode considered is (k x ,k y ) = (tt/4,0). Note that k x>y = 2i:n x ^/L and for this 
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Figure 4: Normalized field intensity of a soft {kt = 108MeV) mode vs proper time t 
(in units of fm) for = 200MeV (squares), /u = lOOMeV (pluses), and /x = 50MeV 
(diamonds). Solid line, nearly coinciding with the = 50MeV curve, is the LPTh 
prediction. 

case, n x = 20, n y = 0. The solid line in the figure is the square of the Bessel function 



Jq(ujt) where uj = y 2(2 — cos(k x ) — cos(k y )). Since the r-direction in our simulation 
is continuous, the time dependence of the high transverse momentum perturbative 
modes should agree with the continuum perturbative result which predicts a time 
dependence proportional to Jq(cjt) for the field intensity of the transverse gauge 
fields. The continuum dispersion relation uj = \k±\ is however modified into the 
lattice dispersion relation shown above. We see from the figure that the anticipated 
agreement between the lattice results and perturbation theory is quite good. 

In Fig. ||, the field intensity of the transverse gauge field normalized by /i 4 at 
r = is plotted as a function of the transverse momentum in physical units. The 
lattice results for the different values of /i described in the caption are compared to 



lattice perturbation theory result-Eq. 69 in appendix B-for the field intensity. The 
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Figure 5: Time history of the energy density in units of fj, 4 for fj, = 200MeV (squares), 
[i = lOOMeV (pluses), and [i = 50MeV (diamonds). Error bars are smaller than the 
plotting symbols. Proper time r is in fm. 

LPTh result (which would be the mini-jet distribution in the continuum) agrees very 
well with the lattice result for small /i upto very small values of kt- However, strong 
coupling effects grow with increasing [i (the lattice size L is fixed here) and we see 
deviations from the perturbative predictions at larger values of k t . This trend is 
enhanced further at larger values of \x than those shown here. The non-perturbative 
effects due to the non-linearities in the Yang-Mills equations seem to temper the 
1/kf behaviour predicted by perturbation theory. Whether this reflects the presence 
of a mass in the theory needs further investigation. 

In Fig. f|, we plot the same quantity as in Fig. 2, but now for three different 
values of /i and for the first non-zero momentum mode (k x , k y ) = (1, 0). The lattice 
size L = 160 is the same as previously. For the smallest value of \x = 0.0177, there 
is again an agreement with the Bessel behaviour predicted by perturbation theory. 
However at the larger values of = 0.035,0.07, one sees significant deviations away 
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from the Bessel behaviour. Indeed, the modes seem to saturate at larger values of 
r. At face value this is unexpected because one expects that as the systems cools at 
late times, even the small k± modes should eventually die out. It may be that we 
have to wait for times much longer than those studied to see this. That would indeed 
be very interesting because these times would be much greater than the natural time 
scale r ~ l/g 2 /J- at which we expect non-linearities to dissipate (see the discussion of 
Fig. H below). 




2 4 6 



Figure 6: Energy per mode k normalized to the total energy as a function of proper 
time in fermis. From top to bottom, the curves correspond to modes k = 0, 108 and 
432 MeV respectively. As in the previous figure, fi = 0.41 GeV. 

A straightforward explanation for this behaviour is that it is some kind of gauge 
artifact. This is because even though the fields satisfy the Coulomb gauge condition 
at t = 0, they no longer do so at later times. Thus one may need to fix Coulomb gauge 
at each "measured" time to properly interpret what we have called as field intensities 
as such. Large modes may be unaffected by the gauge fixing but small modes will 
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be affected. The simplest test of this explanation is to study the correlators of gauge 
invariant quantities at late times. If they display the same "saturation" then the 
effect cannot be dismissed as a gauge artifact. We will return to this point a little 
later. 

Before we do that, we would like to discuss the time dependence of the energy 
density, shown in Fig. [| for different values of \i. At late times, from general con- 
siderations we expect that £ oc 1/t and that is what we see. Indeed, we can see 
qualitatively that the time at which this behaviour is seen is roughly tx 

Thus, even though we clearly see non-perturbative behaviour at low kt, the 
apparent saturation of these modes at large times has no impact on the behaviour 
of the energy density at late times. As fi is decreased, the magnitudes of the energy 
density appear to converge to a fixed value. This behavior of the energy has a 
purely kinematic reason: for any fixed finite rapidity the proper time t asymptotically 
approaches the real time t, with respect to which the energy is conserved. LPTh also 
predicts this behavior. 

We now describe preliminary estimates of the energy distribution, obtained from 
energy-energy correlators, as described in Section 4. The computation was performed 
for an 80 x 80 lattice and for /i = 0.41 GeV. With the small data sample available (50 
independent configurations for the values of proper time r considered), we found that 
the correlation functions C(xt) are consistent with at distances \xt\ > 12 in lattice 
units. We therefore approximated C{xt) as for \xt\ > 12. While the quality of our 
data at this point does not permit a quantitative description of the gluon distribution, 
the data do help us address the question as to whether the unusual behavior of the 
soft modes, as shown in Figure ||, has observable consequences. 

In Fig. ^, we plot the time dependence of the ratio of the energy of the kth mode 
to the total energy of the system (measured directly). The ratio, for the momenta 
considered (0, 108 and 432) MeV, goes to a constant at late times. Since Fig. || clearly 
shows that the energy dies off as 1/r, our result suggests that at late times the energy 
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Figure 7: Energy per mode k normalized to energy of zeroth mode as a function 
of k in GeV units. From bottom to top, the curves correspond to proper times of 
r = 0, 1.4, 2.8, 5.6 fermis. Here y. = 0.41 GeV. 

of the low kt modes must die off as 1/r too. We therefore have an indication that the 
apparent saturation of the field intensities in Fig. ^| is not meaningful and is a gauge 
artifact. 

In Figure |7], we plot the total energy per mode normalized to the energy of 
the zeroth mode as a function of k for the proper times (from bottom to top) r = 
0, 1.4, 2.8, 5.6 fm. Interestingly, the ratio plotted gets flatter at larger values of r. We 
interpret this to mean that the share of the energy in the lower-A^ modes decreases 
at late times. This is another indication of a benign behavior of soft modes at late 
times, contrary to what is suggested by Figure ||. 

6 Summary and outlook 

In this paper, we have performed a non-perturbative study of the production and 
space-time evolution of gluon mini-jets in the central region of nuclear collisions at 
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very high energies. This program was proposed in Ref. p3] and a brief discussion of 



our results can be found in Ref. [55]. 

Our work is based on a classical effective field theory description of the small 
x modes in nuclei at very high energies. This effective theory contains a scale [i 
which is proportional to the large gluon density at small x. The large gluon density 
ensures that even if the coupling is weak, the fields may be highly non-perturbative. 
Since the approach is classical, it is possible to study the real time evolution of 
these non-perturbative modes in a nuclear collision. The approach has the attractive 
feature that it may eventually provide a self-consistent picture of high energy nuclear 
collisions both before and after the nuclei collide. 

At large transverse momenta, kt 3> ctsU, the predictions of the effective theory for 
gluon production should agree with perturbative mini-jet calculations. As discussed 
in Section 5, our lattice results agree with the lattice perturbation theory analogue 
of the continuum mini-jet predictions. At smaller transverse momenta, (in particular 
for large fj,) we see significant deviations from perturbation theory. We also notice 
that the field intensities, \A(kt, t)| 2 , of the small kt modes do not die off at the late 
times studied but appear to saturate. 

It is difficult to interpret the behaviour of the small kt modes since they are 
"off-shell": their dispersion relation u{kt) is non-trivial. Unlike the large momentum 
limit, the field intensities for the small kt modes does not have the simple interpre- 
tation of being the number distribution of produced gluons. As discussed in the 
previous section it is therefore important to look at gauge invariant quantities which 
can be interpreted straightforwardly. In the previous section we discussed a number 
of gauge invariant equal time spatial correlators. These correspond to correlators of 
components of the energy-momentum tensor. In particular, we suggested a gauge 
invariant estimate of the energy distribution on the lattice as a function of time. In 
principle, we should be able to compute the range of energy-energy corelators mea- 
sured in high energy collisions (see Basham et al. |54| and references therein). A more 
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extensive discussion of these correlators will be presented at a later date [47|. 

An important issue not addressed in this work is that of equilibration- do the 
gluons produced equilibrate? We saw in Fig. | that the energy density shows the 
expected behaviour after times oc l/g 2 ^. This however is no proof of thermalization. 
In the classical approach, clearly the high kt modes are not "thermal" . Whether this 
is the case for the softer modes needs to be studied further. It would be interesting 
to relate this approach to that of Miiller and collaborators 1 45 , [49|, ^] . 

We would like to stress that this work is not a quantitative study of gluon pro- 
duction in heavy ion collisions but a qualitative study of non-perturbative effects 
that may be important in the central region of these collisions. Changes that can 
be made in future to bring our results closer to experiment include a) changing the 
gauge group to SU(3), b) relaxing the boost invariance assumption, and c) modifying 
the boundary conditions. Regarding the last point, we have used periodic boundary 
conditions. Modifying these would, for instance, be useful in computing the Poynting 
flux through the nuclear surface. 

Though our choice of gauge A T = is convenient from the point of view of a lattice 
simulation, results in this gauge do not lend themselves easily to a diagrammatic 
interpretation. In this regard it might be useful also to consider a similar computation 
in Coulomb gauge where this is the case [j51{|. Results of such computation will be 
presented at a later date. 
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Appendix A: Numerical method 



In integrating equations of motion of a classical gauge theory it is important to ensure 
that the redundant gauge degrees of freedom do not become observable through 
numerical errors. Hence the integration scheme must respect the Gauss constraints of 
the theory. To this end, the algorithm of choice for the problem at hand is the leapfrog 



algorithm [46]. This integration scheme requires that the Hamiltonian be a sum of 
gauge-invariant kinetic K (dependent on fields only) and potential V (dependent 
on conjugate momenta only) terms, and under this condition, has the advantage of 
respecting the Gauss constraints exactly J45[] . Our lattice Hamiltonian (filD obviously 
has the necessary K + V form, thus the leapfrog algorithm is applicable. 

The light-cone Hamiltonian ([II]) depends explicitly on the proper time r. For 
this reason, the standard leapfrog scheme, designed for time-independent Hamiltonian 
functions, requires some minor adaptation in the present case. In order to explain 
the suitable version of the algorithm, we collectively denote the E, p momentum^ 
variables by P and the U, $ field variables by Q, respectively. The K and V terms 
in the Hamiltonian then have the form K(P,t) and V(Q,t). In the following, the 
argument r of K or V is thought of as a fixed parameter. With this notation, the 
leapfrog step of size A, propagating the system in proper time from r to r + A, has 
the following form. 

1. Integrate the equations Q = {K(P, r), Q} between r and r + A/2. 

3 Strictly speaking, E variables are not conjugate momenta. Indeed, their Poisson brackets with 
each other do not necessarily vanish. However, the scheme as presented only requires that {K, E} = 0, 
and E need not be momentum variables. 
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2. Integrate the equations P = {V(Q, r + A/2), P} between r and r + A. 

3. Integrate the equations Q = {K(P, r + A), Q} between r + A/2 and r + A. 

Here the integration of the equations of motion is assumed to be performed exactly at 
every substep. It can be readily seen that the step as described is exactly reversible, 
i.e., beginning with the final values of P, Q and performing the step with A replaced 
by (—A), one arrives at the initial P, Q. This property of the leapfrog algorithm, 
evident for a time-independent Hamiltonian, is preserved here by suitably choosing 
the explicit proper-time argument of K and V at each substep. The time reversibity 
guarantees that the integration error, obviously of the order no lower than A 2 , is in 
fact 0(A 3 ). 

In order to impose the initial conditions for the proper-time Hamiltonian evolu- 
tion, we first need to solve the lattice Poisson equation (|36|). We do so using the over- 



relaxation method [pq] . The consistency of Eq. 36 requires that the zero-momentum 



component of the color charge density p vanish. We therefore first generate the color 



charge distribution as a normal deviate (35), then subtract from every p q j the spatial 
average J2j Pq,j/ N - 

We fix the Coulomb gauge on initial configurations (see also the discussion in 
section 2.2). The lattice Coulomb gauge reads 



Tr 



E - u'J) - E ~ tfJ-V) ° a 



0, (55) 



where C/j ^ is the link matrix which satisfies the Coulomb gauge condition ( p5[) . We 
use the standard overrelaxation method |37]] for gauge fixing. 

Appendix B: Lattice perturbation theory 

At large transverse momenta, a test of our lattice results is that they agree with 
those of lattice perturbation theory. Here we present a derivation of classical gluon 
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production in lattice perturbation theory. In the continuum limit, we will show that 



our result agrees with Eq. (27) of Ref. [23]. We also obtain an expression for the initial 
kinetic energy on the lattice. We will later compare these analytical expressions to 
the full lattice results at large transverse momenta. 

We restrict our discussion to the gauge group SU(2) and consider the initial 



condition (38) for the link matrix U 

u = (uw + u®)(uw i + uw i y 



Here the labels of the nuclei are written as superscripts. Now recall that since Uj 



1.2 



are pure gauges, they can be written as (see Eq. 34) 



j,n 3 j+n ' 



where V = exp(iAy) and i = 1,2. Hence, 



U® = 1 + 

3,n ' 



i ( A W - A W 

1 V 3 



+ ( a^aW 



3 3+n 2 



(AS J) ) 2 + (A« n ) 2 l)+0(A 3 ). (56) 



The pure gauge solution of the Yang-Mills equations for a single nucleus dictates that 
Vj_AW = (see also Eq. ||). Therefore A ~ 0(fj) and we have kept terms in the 
expansion above up to 0(/j, 2 ). Then substituting Eq. ^ in the expression for U, we 
obtain 



U hn = 1 + %L j>n + - (Q jjn - Qt n - L\ n ) + 



(57) 



where / is the SU(2) identity matrix, 



L 3,n = E ( A . 



W _ A (0 



j+n 



lZ 3,n ~ "j,n ' 



(58) 



i=i 



and 



2 r 



i=l 



A ? A Sn"^((A?) 2 + (4t) 2 



(59) 
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We work in the A T = gauge, supplemented by the Coulomb gauge condition 
( |55| ) at r = 0. The Coulomb gauge-transformed link matrix U' is related to the 
original link U as 

Uj,n' = WjUj, n W} +n . 

This is analogous to the transformation carried out in the continuum perturbation 
theory derivation discussed in section 2.2. For fj, = the Coulomb gauge is satisfied by 
the trivial configuration U = I, Hence the gauge transformation Wj can be expanded 
in powers of small \i: 

Wj = exp (7 + i/jgj + ifi 2 i]j + C(/i 3 )) . 

We now need to determine £ and i] from the Coulomb gauge condition (|55|). Writing 
out U' to order /i 2 , we obtain 

+ CjCj+n + (Lj^j+n - ijhj, n ) + ^(Qj,n ~ <5j >n ~ L "j,n)^ • ( 60 ) 

To lowest order, the Coulomb gauge condition (|60|) implies that 

" 6-+n = — ^ = - (Af + Af ) . (61) 

Now consider the Coulomb gauge condition to order /^ 2 . We need to determine rjj to 
obtain Wj to that order. After some algebra, one can show that 

^41K ),a h +a Sj +(1)h(2) ) ■ (62) 

Above, A is the usual lattice Laplacian 

A(0 = 2 ^ (l-cos(/ n )), 

71=1,2 
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and l n = k n - a. Having computed £ and r] using the Coulomb gauge fixing condition, 
we are now in a position to compute a'— the lattice analog of the field e l in section 2. 
We have 



2ia' hn = Uj, n ' - U jt S = 2i ( Vj - Vj+n ) - ([AW, Ag n ] + (1) <-» (2) J + 0(//) . (63) 

To determine the field intensity we need to compute the lattice Fourier transform 
of the above expression. Using the Fourier transform of A, 

(iV-l)/2 



.(1) A (2) 



A (D,6 



1 



exp(2Trip ■ Xj/L)k\ 



A 



(64) 



P=-(N-l)/2 

after quite some algebra, we obtain the result that 



/a 
a Ln 



c abc 



N 2 



E 

L V 



(1 - e iln ) 



E(-) 2 



l"n' \ ■ ( In' >/ 
sm — sin L/ 

2 J V 2 



(65) 



It is useful to compare this result to the corresponding continuum expression. Here 
and in the following the continuum result is obtained by setting a[ — > a' k a/g, l n — ► kid 
and J2i' ~* L 2 J m^b and letting a — ► 0. In the case of Eq. 65 this prescription gives 



/« _ abc I d 2 k> /k n k n , \r(lU., ,/ \X(2),c 

a fc,n — 6 / J^ZyI \ ~1J2 °V A fc' H'V K n' ) A k-k' 



(66) 



(2tt) 2 V ^ 2 

This result is identical to the Fourier transform of Eq. (27) in Kovner, McLerran and 
Weigert ||]. 

We can now compute the field intensity directly. It is defined as 

Kl = 2^\ a l,n a -I,n/P> 



(67) 



where (• • -) p denotes the Gaussian averaging over the sources. On the lattice, the 
continuum condition VjA® = translates into 



A 



(0 



-\xL 



W 
Vl 

A(J)' 



(68) 
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where rji = p~\a /(JiL. Hence (77, 



Vp 



and {rji-qv 



81 -i>. Substituting Eq. 35 



in Eq. 67 and performing the averaging over p's, we obtain finally the following 
expression for the field intensity 



Kl 2 



3// 4 w A(2f - l)A(l) - [A(V -I)- A(l')f 
1 ^-^ 



(69) 



2A(l)^ A 2 (l>)A 2 (l> -I) 

Here Ylv nieans that terms with V = I or V = are omitted. Taking the continuum 



limit of Eq. |69j, just as we did for Eq. |65| and keeping in mind that \i in this expression 
is given in units of 1/ag 2 (cf. Section 4), one recovers the corresponding continuum 
expression of Ref. [ 23 1 for N c = 2. In the section on our numerical results we compare 



our lattice results for the field intensity to Eq. 69. 

We now compute the initial kinetic energy of the scalar field $ in lattice pertur- 
bation theory. This can also be checked against our lattice results in weak coupling 
and provides yet another test of the numerics. In the continuum, using Eq. ^ and 
Eq. O, we have 



with i = l,2;j = 1, 2. On the lattice, writing 



(70) 



a ■ a 



_Li "_Li 



(71) 



and using a~- \ = (A^(xj + an) — A^(xj))/c 



H (exp(»Z w )- l)m C 2-Kits, 



N 



A(l) 



(72) 



we have (after Gaussian averaging over 77's) 



^ sm(l n ) sm(l n ,) I +ig |^sin 2 (^)sin 2 (^) 



A*(l) 



(73) 
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Let us take a look at the continuum limit of this expression. For n ^ re', the first 
term vanishes. For n = re', this term is oc /x 4 log 2 (L/a) ; it is logarithmically divergent 
in the infrared. The second term is infrared safe for all n,n'. Hence 

p a p a — > A + Blog 2 (L/a), (74) 



where A and B are constants which may be determined from Eq. 73 



Thus far we have not specified precisely what the lattice expansion parameter is. 



In the continuum, perturbation theory is applicable when cts/i/kt <C 1. In Ref. [43 
the lattice expansion parameter for a single nucleus was estimated numerically to be 
g 2 fiL. That this is also the case here can been seen from Eqs. |5(| and |68| . 

Finally, the reader may have noted that our lattice perturbation theory results 
were computed at r = 0. As noted in section 2, in weak coupling the spatial and 
temporal distributions factorize. The spatial distributions of the high momentum 
modes are therefore completely specified at r = 0. 
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